\documentclass[preprint,5p]{elsarticle}
\usepackage{amssymb}
\usepackage{algorithmic}
\usepackage{algorithm}
\newcommand{\theHalgorithm}{\arabic{algorithm}}
\usepackage{subfigure}

\usepackage{amsthm}
\newtheorem{mytheorem}{Theorem}
\newtheorem{mycorollary}{Corollary}
\usepackage{multirow}

\usepackage{url}

\journal{Computer \& Graphics}

\begin{document}

\begin{frontmatter}

\title{Real-Time Ray Casting of Algebraic B-spline Surfaces}

\author{Feifei Wei}
\author{Jieqing Feng\corref{cor1}}
\cortext[cor1]{Correspondence author, jqfeng@cad.zju.edu.cn}
\address{State Key Lab of CAD\&CG, Zhejiang University, 310058, China}
%\author{noname   }
%\author{noname      \corref{cor1}}
%\cortext[cor1]{Correspondence author                       }
%\address{                                                            }

\begin{abstract}
Piecewise algebraic B-spline surfaces (ABS surfaces) are capable of
modeling globally smooth shapes of arbitrary topology. These can be
potentially applied in geometric modeling, scientific visualization,
computer animation and mathematical illustration. However, real-time
ray casting the surface is still an obstacle for interactive
applications, due to the large amount of numerical root-findings of
nonlinear polynomial systems that are required. In this paper, we
present a GPU-based real-time ray casting method for ABS surfaces.
To explore the powerful parallel computing capacity of contemporary
GPUs, we adopt iterative numerical root-finding algorithms, {\it
e.g.}, Newton-Raphson and regula falsi algorithms, rather than
recursive ones. To facilitate convergence of Newton-Raphson or
regula falsi algorithm, their initial guesses are determined through
rasterization of the isotopic isosurface, and the isosurface is
generated based on regular criteria for surface domain subdivision.
Meanwhile, polar surfaces are adopted to identify single roots or to
isolate different roots, {\it i.e.}, ray and surface intersections.
As an important geometric feature, the silhouette curve is
elaborately computed to floating point accuracy, which can be
applied in further anti-aliasing processes. The experimental results
show that the proposed method can render thousands of piecewise
algebraic surface patches of degrees 6-9 in real time.
\end{abstract}

\begin{keyword}
Algebraic B-spline Surfaces \sep Newton-Raphson method
\sep regula falsi \sep polygonization \sep polar surface \sep
silhouette
%% keywords here, in the form: keyword \sep keyword
\end{keyword}

\end{frontmatter}

%%
%% Start line numbering here if you want
%%
% \linenumbers

%% main text
\section{Introduction}\label{sec:Intro}
Because of its advantages in intersection, blending, offset, and
arbitrary topology, the algebraic polynomial representation is an
alternative to parametric surfaces in geometric modeling. The
piecewise algebraic surface in terms of B-splines, which is
abbreviated as ABS in this paper, was first introduced by
Patrikalakis and Kriezis~\cite{Patrikalakis89}. The ABS surface is
capable of modeling globally smooth shapes of arbitrary topology
without explicitly specifying smoothing constraints as in
A-patch~\cite{Bajaj95} or NURBS surface~\cite{Che05} modeling.
Compared to subdivision surfaces, which can also model globally
smooth shapes of arbitrary topology, the ABS surface has an
analytical expression, which is important to perform geometry
computations and shape interrogation in geometric modeling, {\it
e.g.}, evaluations of normal, curvatures, offset, intersection.
J\"{u}ttler and Felis~\cite{Juttler02}, and Tong {\it et
al.}~\cite{Tong06} present ABS surface fitting algorithms for
reconstruction of scatter data. However, shape modeling via ABS
surfaces is beyond the scope of this paper.
%The ABS surfaces is an extension of piecewise algebraic surface
%patch~\cite{Sederberg85} in terms of Bernstein polynomials defined
%in a tetrahedron, also called A-patch.

Unlike its parametric counterpart, high quality ABS surface display
is not trivial. The ray casting approach can reveal elaborated
surface topology and geometric details, thus, it is free of sampling
artifacts in pixel accuracy. Furthermore, only polynomial
coefficients of the surface are required to use a ray casting
algorithm, which leads to lower memory and bus bandwidth consumption
than does the use of prevalent triangle rendering pipelines.
However, the ray casting approach involves a large number of ray and
surface intersections. As far as we know, there are no real-time ray
casting methods for ABS surfaces.

In this paper, we present a GPU-based real-time ray casting method
for ABS surfaces via iterative root-finding algorithms, {\it e.g.},
Newton-Raphson and regula falsi. The initial guess of the
root-finding algorithm is first given via rasterization of an
isosurface obtained from an isotopic ABS surface polygonization. The
first root is isolated by analyzing configurations between ABS
surfaces and their polar surfaces. Finally, Newton-Raphson algorithm
or regula falsi algorithm is used to compute the first hit point.
Meanwhile, the polar surface is also employed to compute silhouette
points to floating point accuracy. And the silhouette points can be
applied for further anti-aliasing. The hit point is computed in
object space, and thus the overhead of functional composition
between the ray and surface is avoided. Furthermore, the surface
normal for shading is a byproduct of the iterative root-finding
algorithm. The proposed method can display an ABS surface of degree
6-9, which is composed of thousands of piecewise patches, in real
time.
%the polar surface is also employed to compute silhouette points to floating point accuracy,and the silhouette points are used for further anti-aliasing.


The contributions of the proposed method are summarized as follows:
\begin{itemize}
% \item Real-time rendering of piecewise ABS surface of degree above 4 on GPU;
  \item Real-time ray casting of ABS surfaces on GPU;
  \item Rather than solving a univariate polynomial equation, we iteratively solve a nonlinear tri-variate
  polynomial system in object space to find the first hit point, where an initial
  guess is determined via isotopic polygonization;
%  \item Adopting initial guesses by rasterization of isotopic polygonal mesh;
  \item Roots isolation by using polar surfaces;
  \item Explicitly computing the silhouette and anti-aliasing along the silhouette;
\end{itemize}

In the following, we first review related work in Section
\ref{sec:related}. Then we describe the ABS surface in Section
\ref{sec:ABS} and our algorithm in Section~\ref{sec:alg}. Isotopic
polygonization and applications of the polar surface are described
in Section~\ref{sec:MC} ,~\ref{sec:Silhouette} and
~\ref{sec:isolate}, respectively. Details of our implementation are
given in Section~\ref{sec:Imple}. Numerical examples, comparison and
discussion are given in~\ref{sec:result}. Finally we draw the
conclusion and propose future work.


\begin{figure*}[ht]
  \centering
  \includegraphics [width=0.9\textwidth]{flow2}
  \caption{Flow chart of real-time rendering of ABS surfaces}\label{fig:flow}
\end{figure*}

\section{Related work}\label{sec:related}

Currently, there are two kinds of methods used to display an
algebraic surface: sampling it and displaying proxy geometries, {\it
e.g.}, triangles or particles and directly ray tracing or ray
casting it. Displaying an algebraic surface via discrete proxy
geometries is more suitable for prevalent triangle rendering
pipelines. With the rapid development of graphics hardware,
parallelism of triangle rendering pipelines is exploited thoroughly.
Even for very high throughput, the state of the art rendering
pipelines are able to produce stunningly-realistic synthesis images
in real time. The prevalent approaches to approximate or sample an
algebraic surface are Wyvill's polygonization \cite{Wyvill86},
Marching Cubes \cite{Lorensen87}, Bloomenthal's
polygonization~\cite{Bloomenthal94} and the particle system
approach~\cite{Witkin94}. However, sampling artifacts are the main
flaw of these approaches. Level-of-detail techniques can be employed
to alleviate the artifacts to some extent, which will cause
additional overhead. A comparison between ray casting and
polygonization is shown in Fig.\ref{fig:polygonization}. The
polygonization resolution of the left surface is up to $200^3$,
however, it is still plagued by sampling artifacts.

\begin{figure}[h]
\centering
\includegraphics[width=0.45\textwidth]{poly}
\caption{Comparison between our methods and polygonization of the
surface in Fig.\ref{fig:SingularSharp}. The yellow surfaces are
rendered by polygonization, while green by our method.The examples
illustrate that the ray casting images can always reveal the details
of the ABS surface more elaborately, such as the singular point and
the sharp edge, while the polygonization may miss these fine
features.} \label{fig:polygonization}
\end{figure}

Although ray casting can address the above problems, it usually
cannot meet the demands of real-time applications. Its most
time-consuming procedure is to find the first hit point between a
ray and the surface. In general, this is a root finding problem of
polynomials. Analytical solutions exist for polynomials of degree 4
or lower. According to this observation and some fundamental work
in~\cite{Blinn05}, Loop and Blinn \cite{Loop06} proposed a real-time
ray casting algorithm for A-patches up to degree 4. However, it is
not trivial to apply this approach to render ABS surfaces directly
since the configurations of a cube projection on the screen plane
are more complex than those of a tetrahedron. It is also not
efficient to process overlapped domain cubes with complex depth
relationships.

In the past decades, many numerical root-finding algorithms have
been employed to ray cast parametric or algebraic surfaces of
degrees above 4. Kajiya~\cite{Kajiya82} reduced ray tracing of
parametric surfaces to root finding of a high degree polynomial
using Laguerre's method. Hanrahan~\cite{Hanrahan83} implemented ray
tracing of algebraic surfaces by using Descartes rule of signs for
root isolation and Newton's method for root refinement. Kalra and
Barr \cite{Kalra89} employed a Lipschitz constant and surface
gradient to guarantee correct ray and surface intersection. Interval
arithmetic constructs a convex hull around a function over a given
domain. By using this property, Mitchell~\cite{Mitchell90} isolated
the roots using repeated interval bisection until an interval
contains a single root. B\'{e}zier clipping, proposed by Nishita
{\it et al.}~\cite{Nishita90}, finds the roots of a Bernstein
polynomial equation based on the convex hull property. Recently,
M{\o}rken and Reimers~\cite{Mrken07} presented an unconditionally
convergent root finding method via knot insertion for splines.

Utilizing the GPU, it is possible to achieve real-time ray tracing
of algebraic surfaces. Pabst{\it et al.}\cite{Pabst06} and Kanamori
{\it et al.}\cite{Kanamori08} implemented an iterative version of
B\'{e}zier clipping with a fixed-size stack on GPU for ray casting
trimmed NURBS surfaces and metaballs respectively. Knoll {\it et
al.} proposed an interval bisection algorithm for real-time ray
casting general algebraic surfaces via SSE instruction-level
optimization~\cite{Knoll07} and GPU~\cite{Knoll09} respectively.
Both interval arithmetic and B\'{e}zier clipping require recursive
subdivision until a tolerance is satisfied. Interval evaluation of
Bernstein polynomials is numerically instable due to large numbers
of MAD (addition and multiplication) operations on GPU. Meanwhile,
Rounded Interval Arithmetic is conservative~\cite{Maekawa93}, but is
inefficient on GPU. Reimers and Seland \cite{Reimers08} performed a
functional composition between ray and surface via polynomial
interpolation and employ B-spline knot insertion for root-finding.
The root-finding algorithm is iterative and more suitable for the
SIMD architecture of GPU than is the B{\'e}zier clipping approach.
However, functional composition is an additional overhead. Singh and
Narayanan \cite{Singh10} proposed an adaptive marching point method
to find the first hit point, wherein the step size is determined by
distance from the surface and its closeness to the silhouette.
Although it is very efficient on GPU, root isolation via rule of
signs is not always robust. Furthermore threshold estimations of
step size are not a trivial work for piecewise surfaces, since they
are patch dependent.

\section{Algebraic B-spline surfaces}\label{sec:ABS}

Let $\mathbf{X}=[x_{0},x_{1},\cdots,x_{M+m+1}]$,
$\mathbf{Y}=[y_{0},y_{1},\cdots, y_{N+n+1}]$ and
$\mathbf{Z}=[z_{0},z_{1},\cdots,z_{Q+q+1}]$ be three non-decreasing
knot vectors along the $x$, $y$ and $z$ directions. An ABS surface
is defined as:

\begin{equation}\label{equ:ABS}
\hspace{-8mm}
F_{ABS}(x,y,z)=\sum_{i=0}^{M}\sum_{j=0}^{N}\sum_{k=0}^{Q}w_{ijk}N_{i}^{m}(x)N_{j}^{n}(y)N_{k}^{q}(z)=0
\end{equation}
where $N_{i}^{m}(x),~N_{j}^{n}(y)$ and $N_{k}^{q}(z)$ are the
B-spline basis functions of degree $m$, $n$ and $q$, defined by the
knot vectors $\mathbf{X}$, $\mathbf{Y}$ and $\mathbf{Z}$. The
scalars $\{w_{ijk}\}$ are weight coefficients. They have similar
functions to control points in their parametric counterpart and can
be used to adjust the shape of the ABS surface. In this paper, the
degree of the ABS surface along each direction is no greater than 3.
They can model the $C^0$, $C^1$ or $C^2$ smooth shapes, which means
that the shape has common position, tangent and curvature values
along the boundaries of the patches respectively. In general, they
can meet the requirements in the practical applications, {\it e.g.},
$C^0$ or $C^1$ smooth shapes for computer animation and computer
games, and $C^1$ or $C^2$ shapes for automotive design. By using the
knot insertion algorithm of B-splines, the ABS surface defined on a
knot box $[x_{i},x_{i+1}]\times[y_{j},y_{j+1}]\times[z_{k},z_{k+1}]$
$(m\le i \le M, n\le j \le N, q\le k \le Q)$ can be converted into
an algebraic surface patch in terms of tensor product Bernstein
polynomials:
\begin{equation}\label{equ:ABez}
\hspace{-5mm}
\begin{array}{l}
  F(x,y,z)= \\
  \displaystyle\sum_{i=0}^{m}\sum_{j=0}^{n}\sum_{k=0}^{q}f_{ijk}B_{i}^{m}(u(x))B_{j}^{n}(v(y))B_{k}^{q}(w(z))=0\\
\end{array}
\end{equation}
where
%\begin{equation}\label{equ:3}
\[
\hspace{-5mm}
u(x)=\frac{x-x_{i}}{x_{i+1}-x_{i}},~v(y)=\frac{y-y_{j}}{y_{j+1}-y_{j}},~w(z)=\frac{z-z_{k}}{z_{k+1}-z_{k}}.
\]
%\end{equation}
$f_{ijk}$ are coefficients similar to $w_{ijk}$ of ABS surfaces, and
$B_{i}^{m}(u)$, $B_{j}^{n}(v)$ and $B_{k}^{q}(w)$ are Bernstein
basis functions of degree $m$, $n$ and $q$. As a result, we can
obtain $(M-m+1)\times(N-n+1)\times(Q-q+1)$ piecewise algebraic
surfaces in Bernstein form.

\section{Rendering algorithm}\label{sec:alg}

\subsection{Algorithm Overview}\label{sec:overview}

In general, ray and surface intersection can be reduced to a
root-finding problem of a univariate polynomial equation via
functional composition. Functional composition between an ABS
surface and a ray is a complex procedure \cite{DeRose93}, and is
impractical for the SIMD architecture. In our setting, ray and
surface intersection is performed in object space as a root-finding
problem of a tri-variate nonlinear system as follows:
\begin{equation}\label{equ:RaySr}
%\[
\left\{ {\begin{array}{ccc}
   %{F = \sum\limits_{i = 0}^M {\sum\limits_{j = 0}^N {\sum\limits_{k = 0}^Q {w_{ijk} N_{i}^{m}(x)N_{j}^{n}(y)N_{k}^{q}(z) = 0} } } }  \\
   F(x,y,z) &=& 0  \\
   a_h x + b_h y + c_h z + d_h  &=& 0  \\
   a_v x + b_v y + c_v z + d_v  &=& 0  \\
\end{array}} \right.
%\]
\end{equation}
where $F$ is a piecewise tensor-product Bernstein polynomial
described in Equation~(\ref{equ:ABez}), and the two linear equations
are the horizontal and vertical scan planes that pass both the
viewpoint and current pixel, respectively. Reimers and Seland
adopted functional composition to reduce the ray and the algebraic
surface intersection to a root-finding problem of univariate
polynomial \cite{Reimers08}, which may contain many evaluations of
the algebraic surface. In our approach, Equation~(\ref{equ:RaySr})
is solved in the object space directly. Thus additional
computational costs for functional compositions are saved.

%In this way, the functional composition between the ray and the
%function, which contains many evaluations of the ABS surface, is
%avoided.

The nonlinear system (\ref{equ:RaySr}) can be solved via
Newton-Raphson method, which can converge quadratically if a good
initial guess is given. A rough initial guess can be obtained via
rasterization of an isosurface of the underlying ABS surface, which
is an isotopic polygonization as stated in Section~\ref{sec:MC}.
However, the initial guess cannot guarantee convergence of the root
finding algorithm near the silhouette due to ill-conditioned
computations. To overcome this problem, the silhouette, which is the
intersection between the underlying ABS surface and its polar
surface, is processed specially and robustly. For regions near the
silhouette, the polar surface is employed again to isolate an
interval containing the first hit point. The regula falsi method is
adopted to find the first hit point in this interval. For other
regions, the Newton-Raphson method can achieve quadratic convergence
robustly. The flow chart for the proposed method is shown in
Fig.\ref{fig:flow}, and its description is given in
Algorithm~\ref{alg:flow}.


%\begin{enumerate}\label{algorithmstep}
%  \item Extract isosurface with Marching Cubes, where the cube resolution satisfies that the isosurface mesh is a topological equivalent to the original
%ABS surface;
%  \item Render the isosurface in GPU pipeline, and obtain the initial
%  estimation of Newton-Raphson method;
%  \item Compute the silhouette via polar surface, and rectify the initial estimation nearby the silhouette;
%  \item Sort the initial estimation based on patch index, conduct Newton-Raphson iteration on the initial estimation cluster;
%  \item Refine the boundary where adjacent patches meet, compute the shading result and anti-aliasing.
%\end{enumerate}




\begin{algorithm}[ht] \caption{}
\begin{algorithmic}[1]\label{alg:flow}

\STATE {\textbf{Input:} ABS surface}

\STATE {Isotopic polygonization of ABS surface;}

\STATE {For each patch of ABS surface, check existence of its polar
surface via signs of control coefficients;}

\REPEAT

\STATE {Peel \textit{i}-th layer of visible isosurface rendering
result as initial guess;}

\STATE {Compute silhouette points and rectify the initial guess near
the silhouette;}

\STATE {\textbf{Scan} the initial guess and \textbf{Sort} them by
patch index to generate pixel clusters;}

\FOR {each pixel cluster}

\IF {there is no polar surface}

\STATE {Compute first hit point between ray and patch via
Newton-Raphson method by using the initial guess;}

\STATE {Record surface normal in last iteration;}

\ELSE

\STATE {Compute first hit point between ray and polar surface by
using Newton-Raphson method;}

\STATE {Compute first hit point between ray and patch via regula
falsi method in the interval bracketed by ray entry and polar
surface;}

\STATE {Calculate surface normal;}

\ENDIF

\ENDFOR

\UNTIL {No intersection is available.}

\STATE {Shading and anti-aliasing.}

\end{algorithmic}
\end{algorithm}



\subsection{Isotopic Polygonization}\label{sec:MC}

Topologically consistent polygonization of an implicit surface is
not trivial. Several algorithms have been proposed to address this
problem, which are based on Morse theory, stratification theory,
{\it etc}\cite{Barton97, Alberti05, Alberti09}. Because the ABS is
transformed to Bernstein basis representations in rectangular grids
as in Equation (\ref{equ:ABez}), and Alberti {\it et.al.}
\cite{Alberti05}proposed a regular criteria according to a
generalization of Descartes' rule that can guarantee that the
resulting mesh is isotopic to the algebraic surfaces. We keep
subdividing the representation until the surface in each box boils
down to the case where the implicit object is isotopic to its linear
approximation in the cell. Thus, we can estimate the maximum
subdivision resolution of the ABS surfaces conservatively by
adopting the maximum resolution among all the algebraic patches.
Then the ABS surface is polygonized in parallel on the GPU according
to this resolution. Due to the tensor-product property of the ABS
surface and its low degree (less than 4 along each direction), we
augment the polygonization by use of accurate vertices in the
isosurface. The accurate vertex, {\it i.e.}, the intersection
between the ABS surface and voxel edge, can be analytically
evaluated by using Blinn's efficient GPU algorithm~\cite{Blinn05}.
The additional cost of augmented polygonization is less than 10ms in
general, according to our implementations. However, accuracy
improvement in polygonization will be helpful to enhance the
performance of ray and surface intersection computation.

As indicated by Alberti {\it et al.}~\cite{Alberti05}, the above
isotopic polygonization is confined to a smooth surface with finite
singular points, where a singular point is defined
as$f=f_x=f_y=f_z=0$. Thus we also assume that the ABS surface has
finite singular points in this paper. The singular points can be
found robustly via the Projected Polyhedron
algorithm~\cite{Patri02}. The isosurface near a singular point is
described as a cone-like shape~\cite{Alberti05}. We suppose that the
ABS surfaces are potentially more suitable to model smooth shapes,
rather than shapes with sharp features such as infinite singular
points or self-intersections.

For the same viewpoint as that of the following recasting, the
isosurface is first sent to the traditional polygonal rendering
pipelines. Then, a rasterized image together with its z-buffer is
obtained. Each z-value in the z-buffer together with its
corresponding pixel $(x, y)$ is coalesced as a coarse initial guess
for the following root-finding iteration in the raycasting.

%We render the isosurface at the same viewpoint via triangle
%rendering pipelines. Each $z$-value in the z-buffer together with
%its corresponding pixel $(x,y)$ are coalesced as a coarse initial
%guess for the following root-finding iteration.





\subsection{Silhouette and Polar Surfaces}\label{sec:Silhouette}
The silhouette plays an important role in the high quality display
of an algebraic surface. Sederberg and Zundel~\cite{Sederberg89}
first proposed the concept of the polar surface to compute the
silhouette. Its definition is described in homogenous coordinate
space. For an algebraic surface $F(x,y,z,w)=0$ of degree $n$ in
homogeneous coordinate space, its polar surface $P(x,y,z,w)$ with
respect to viewpoint $E(x_e,y_e,z_e,w_e)$ is of degree $n-1$ and is
defined as:
\begin{equation}\label{equ:Polar}
    P(x,y,z,w)=x_{e}F_{x}+y_{e}F_{y}+z_{e}F_{z}+w_{e}F_{w}=0.
\end{equation}


At this point the polar surface is rewritten in Euclidean space by
setting $w=w_e$. If the projection is orthogonal, $w_e=0$; if
perspective, $w_e=1$. For the sake of convenience, the polar surface
is still noted as $P(x,y,z)=0$. The degree of the ABS surface in
Equation~(\ref{equ:ABS}) is $(m+n+q-1)$. To facilitate the following
computations on GPU, its degree is elevated to $m\times n\times q$.
The silhouette curves with respect to a viewpoint is defined as the
intersection between the ABS surface and its polar
surface~\cite{Sederberg89}. We then perform scan conversion of
silhouette points horizontally or vertically according to the
following system:
\begin{equation}\label{equ:RaySrPo}
\left\{
\begin{array}{c}
    F(x,y,z)=0\\
    P(x,y,z)=0\\
    a_{k}x+b_{k}y+c_{k}z+d_{k}=0,k=v~or~h
\end{array}
\right.
\end{equation}
where the third linear equation represents the scan plane.

%\begin{figure}
%\centering\subfigure[]{
%\includegraphics[width=0.44\textwidth]{Classification.eps}
%\label{fig:Classification}}
%\subfigure[]{
%\includegraphics[width=0.44\textwidth]{edge.eps}
%\label{fig:Edge} }\caption{Computation of silhouette.
%\subref{fig:Classification} Classifications of silhouette points,
%\subref{fig:Edge} Padding of the gaps nearby the silhouette.}
%\end{figure*}

\begin{figure}[h]
\centering
\includegraphics[width=0.4\textwidth]{Classification.eps}
\caption{Classifications and iteration of silhouette computation.
The horizontal and vertical arrows indicate the horizontal and
vertical scan planes. The initial guesses in the gaps near the
silhouette are obtained via linear
interpolation.}\label{fig:Classification}
\end{figure}

Equation~(\ref{equ:RaySrPo}) is solved via the Newton-Raphson method
on GPU as indicated in the sixth step of Algorithm~\ref{alg:flow}.
The silhouette of the isosurface can be adopted as its initial
guess, namely the silhouette initial guess, via edge detection of
its rasterized image. To perform a well-conditioned iteration, the
horizontal or vertical scan plane should be selected carefully in
Equation~(\ref{equ:RaySrPo}). The gradient vector of $F(x,y,z)$ at
the silhouette initial guess is projected onto the screen plane,
which is noted as $(n_x, n_y)$. If $ |n_x| \approx |n_y|$, both the
horizontal plane and vertical scan plane should be selected as shown
in Fig.\ref{fig:Classification}. There are two silhouette points
obtained, which is important for the convex silhouette. If $ |n_x|
> |n_y|$, the vertical scan plane should be selected; otherwise, the
horizontal scan plane should be chosen. If a silhouette curve
segment is convex, then there will be a gap region between the
isosurface silhouette and surface silhouette. In this gap region,
there is no initial guess definition for ray and surface
intersection. The missed initial guess can be obtained via linear
interpolation of the initial guess at the isosurface silhouette and
the actual silhouette point.

\subsection{Root isolation via polar surface}\label{sec:isolate}

To compute the intersection between a ray and the ABS surface, we
solve Equation~(\ref{equ:RaySr}) via Newton-Raphson method by using
the initial guess, which is obtained from isosurface rasterization
in Section~\ref{sec:MC} or linear interpolation in
Section~\ref{sec:Silhouette}. However, near the silhouette, a ray
and the ABS surface tends to intersect at two points $P_1$ and
$P_2$, as shown in Fig.\ref{fig:isolate}. Thus Newton-Raphson
iterations with the initial guess $P_m$ may not converge to the
first hit point $P_1$. To address this problem, we first isolate the
first hit point by using the polar surface. Then, the regular falsi
method is employed to compute the first hit point. The details of
root isolation via polar surfaces are provided in Appendix. Here we
give two corollaries of the polar surface:

\begin{mycorollary}\label{corollary:1}
If the polar surface $\{(x,y,z):P(x,y,z)=0\}$ is null, the ray will
intersect the surface at one point at most, {\it i.e.},
Equation~(\ref{equ:RaySr}) has at most one root.
\end{mycorollary}

\begin{mycorollary}\label{corollary:2}
The intersection point between the ray and the polar surface can
isolate two adjacent intersection points between the same ray and
the corresponding algebraic surface.
\end{mycorollary}

\begin{figure}[h]
\centering
\includegraphics[width=0.40\textwidth]{derivate}
\caption{Root isolation using the polar surface near the
silhouette.}\label{fig:isolate}
\end{figure}

As stated in Algorithm~\ref{alg:flow}, two iterative root-finding
methods are adopted: the Newton-Raphson method and a two-phases
method (root isolation plus regula falsi). According to
Corollary~\ref{corollary:1}, if there is no polar surface defined in
${\mathbf \Re}^3$, a ray and the ABS surface will intersect in one
point at most. In this case, the Newton-Raphson method with a rough
initial guess will converge quadratically. Otherwise, the initial
guess will be close to the silhouette as shown in
Fig.\ref{fig:isolate}, and the Newton-Raphson iterations will be
ill-conditioned due to a near singularity of the Jacobian matrix of
Equation (\ref{equ:RaySr}). To overcome this problem, we isolate the
first hit point $P_1$ of the ray and the ABS surface by using the
polar surface according to Corollary~\ref{corollary:2}. Because the
ray tends to be orthogonal to the polar surface near the silhouette,
we can compute the ray and the polar surface intersection $p_p$ by
using the Newton-Raphson method. Let the ray and the patch domain
intersect at two points, {\it i.e.}, $p_{entry}$ and $p_{exit}$. The
interval $[p_{entry}, p_p]$ or $[p_p, p_{exit}]$ is selected as the
initial interval for the regula falsi method, in which $F(x,y,z)$
has opposite signs at the interval ends. The regula falsi method is
robust, as the first root has been bracketed, even for the pixels
near a singular point. An example is shown in Fig.\ref{fig:surPol}.


\begin{figure}[h]
\centering
 \subfigure[front view]{
\includegraphics[width=0.22\textwidth]{Polar2.eps}\label{fig:subPolar1}
}\subfigure[side view]{
\includegraphics[width=0.22\textwidth]{t_05.eps}\label{fig:subPolar2}
} \caption{An hourglass shape surface with one singular point and
its polar surface. The polar surface isolates two intersection
points of the ray and surface.}\label{fig:surPol}
\end{figure}


\subsection{Shading and Anti-aliasing}\label{sec:Anti}

The surface normal is important for illumination calculations. In
functional-composition-based ray casting approaches, the surface
normal is computed via an additional procedure, either analytically
or approximately. The root finding procedure in our approach is
performed in object space thoroughly, and the extra overhead of
normal computation can be avoided. For example, in a Newton-Raphson
iteration, the surface normal $\nabla F$ is computed in each step.
The $\nabla F$ computed in the final iteration can be adopted as the
surface normal.

In Section~\ref{sec:Silhouette}, the silhouette points have been
computed to floating point accuracy. In screen space, the slope at
each silhouette point can be obtained via projection of the normal.
In this way, each silhouetting pixel can be anti-aliased by using a
percentage coverage filtering methodology.


\section{Implementation on GPU}\label{sec:Imple}
The proposed algorithm is designed for GPU with SIMD architecture.
Except for isosurface rendering via a traditional graphics pipeline,
the proposed algorithm is implemented completely in NVIDIA CUDA.
CUDA provides a C language interface for general-purpose programming
on the GPU. It exposes many important parallel computing
capabilities of the GPU so as to facilitate data-parallel
computations, {\it e.g.}, an efficient scatter-gather mechanism. The
common operations, such as scan, reduce and sort, have been released
in CUDPP\cite{Harris07}, which can be reused conveniently as shown
in Algorithm~\ref{alg:flow}.

%By utilizing the parallel power of GPU,
%the isosurface extracting with proper grid resolution under the Lipschitz condition can be executed in real time.
%So, the performance will not be greatly influenced for dynamical surfaces.

The major difference between the rendering of single algebraic
surfaces and ABS surfaces is that ABS surfaces have a large number
of control coefficients. The manner in which the piecewise surface
patch data is accessed seriously affects the performance of the
rendering algorithm.

According to the specification of CUDA, allocating more threads per
block is better for efficient time slicing. However, the more
threads per block, the fewer registers are available per thread due
to limited register resources on GPU. As a result, the coefficients
of each surface patch and its polar surface could cause high memory
throughput. If we implement the Newton-Raphson method trivially, it
will load more coefficients when the iterations jump from one patch
onto its neighbor. Thus it is impractical to load surface patch data
into registers for every pixel. On the other hand, the simultaneous
memory access of adjacent pixels could not be coalesced into a
single memory transaction, which would result in a performance
penalty.

To address this problem, we propose a patch-index based clustering
scheme. First, each initial guess is associated with a patch index
according to its position. We sort and cluster the initial guesses
by their associated patch index. Then we launch these clusters in
blocks. In this way, threads in each block can share the same patch
data, which can be loaded into the on-chip shared memory. This
scheme provides efficient memory access as fast as registers,
meanwhile limited registers are saved and could be applied to deeper
parallelism.

During the implementation of the Newton-Raphson or regula falsi
method, the iteration could cross the boundary between two patches.
In this case, we still use current patch data via surface patch
extrapolation. After the iteration terminates, we scan the pixels
and cluster those that actually converge to the adjacent patch. Then
we refine these with additional iterations. Such pixels are
typically found near the patch boundary, and thus the pixel number
is limited in general. The additional iterations can also be
performed efficiently even if we read the patch data through texture
memory fetching. The operations of ``scan'' and ``sort'' are both
implemented by using the CUDPP library on GPU.

\section{Results and Comparison}\label{sec:result}
\begin{figure}[h]
     \centering
    \subfigure[1 patch, $G^0$, 57.62 fps]
    {
        \includegraphics[width=0.25\textwidth]
        {drop}
        \label{fig:drop}
    }
    \subfigure[3 patches, $G^0$, 97.65 fps]
    {
        \includegraphics[width=0.18\textwidth]
        {stamp}
        \label{fig:stamp}
    }
    \caption{(a) is an ABS patch of degree 7 with a singular point constructed from the ``Dingdong" surface.
    The sharp edge of (b) is constructed via multiple knots.}
    \label{fig:SingularSharp}
\end{figure}


The algorithm has been implemented on a PC with an Intel Core2 Q9550
2.83GHz CPU and an NVIDIA GeForce GTX 570 GPU. The rendering
resolution is $1024\times768$, and the ABS surfaces roughly fill the
screen. The termination criteria in the root-finding of
Equation~(\ref{equ:RaySr})is that the value of $|F(x,y,z)|$ or the
iterative step of consecutive iterations is less than
$1.0\times10^{-6}$, which is the threshold of single precision
floating point arithmetic.

\begin{figure*}\label{fig:5shapeMorf}
\centering
\includegraphics[width=0.95\textwidth]{5shapeMorf}
\caption{Morphing between two ABS surfaces of degree 6 with
different topologies. The morphing is generated by interpolating the
control coefficients of the source and destination surfaces on the
fly. This can be achieved at approximately 29-35 fps with a
$1024\times768$ resolution.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`}

%(The demo video can be found at
%{http://www.cad.zju.edu.cn/home/jqfeng/RealtimeRaycastingABS.htm}).}
\end{figure*}

We constructed several examples using the ABS surfaces. The shapes
in Fig.\ref{fig:surPol} and \ref{fig:drop} both have one singular
point, while Fig.\ref{fig:stamp} has sharp edges. Fig.\ref{fig:333}
is $C^2$ smooth, and all the others are $C^1$ smooth in all
directions. Their topologies vary, {\it e.g.}, genus 0 in
Fig.\ref{fig:bishop}, genus 1 in Fig.\ref{fig:04} and
\ref{fig:torus}, genus 5 in Fig.\ref{fig:333} and even genus 12 in
Fig.\ref{fig:05}. Some of these are regular shapes, as the torus and
revolution surface in Fig.\ref{fig:L06}; some are free-form shapes,
{\it e.g.}, bunny and knot in Fig.\ref{fig:demo3}. From these
examples, we conclude that the ABS surface can not only represent
globally-smooth shapes, but is also able to model shapes with sharp
features. However, the modeling of the ABS surfaces goes beyond the
scope of the paper and is omitted here.



%Figure~\ref{fig:demo} and Figure \ref{fig:demo3} show some globally
%smooth ABS surfaces of degree 6-9, where different colors indicate
%different patches in Figure \ref{fig:demo}. And we can also create
%ABS surfaces from the signed distance fields of polygonal meshes, as
%shown in Figure \ref{fig:demo3}. It demonstrates that ABS surface is
%capable of modeling globally smooth shape of arbitrary topology.

Because the algebraic surface is capable of processing a dynamic
topology, an example of topological morphing between two ABS
surfaces is given in Fig.\ref{fig:5shapeMorf}, where the source and
the destination ABS surfaces have genus $0$ and $5$, respectively.
Because the proposed algorithm is purely on-the-fly, it can be
applied to render dynamic ABS surfaces and explore the topological
changes.

Information regarding these ABS surfaces, the subdivision resolution
of the isotopic polygonization, running time statistics and
comparisons are listed in Table~\ref{table:patch}. In our algorithm,
only the maximum subdivision resolution and singular point detection
in the isotopic polygonization are evaluated by the CPU. The
successive polygonization is computed once or on-the-fly for the
static or dynamic ABS surface on GPU respectively. All remaining
parts of the algorithm are performed on-the-fly on GPU.

%The numbers of the patches are determined by the signs of control
%coefficients. In the ``Polygonization'' column, ``sub'' means the
%subdivision resolution of isotopic polygonization, ``accurate'' and
%``lerp'' indicate the polygonization time using accurate
%root-finding and linear interpolation to find the edge points.



\begin{table*}[ht]
\centering \caption{The ABS surfaces information and a comparison of
different methods, where ``segments'' is the number of valid domains
of the ABS surface when transformed into Bernstein form as in
Equation (\ref{equ:ABez}), {\it i.e.}, $(M-m+1)*(N-n+1)*(Q-q+1)$,
``$\#$(patches)'' is the number of algebraic patches in terms of
Bernstein basis, ``sub'' refers to the subdivision resolution of the
isotopic polygonization, ``acc(ms)'' and ``lerp(ms)'' represent the
polygonization run time where the vertices of the mesh are computed
Accurately- or Linearly-approximated respectively.}
\label{table:patch}
%\small\addtolength{\tabcolsep}{-5pt}
\begin{tabular}{cccc|ccc|ccc|cc}

  \hline
%  \multirow{3}{*}{ABS}& \multirow{3}{*}{degree}& \multirow{3}{*}{knot segments}& \multirow{3}{*}{patches} & \multirow{3}{*}{sub} & \multicolumn{2}{c}{Marching Cubes(ms)}& \multicolumn{2}{c}{Frame rate(fps)} \\
%  \cline{6-9}
%  &&&&& accurate & lerp  & Ours & RW\\
  \multirow{2}{*}{Fig}& \multirow{2}{*}{degree}& \multirow{2}{*}{segments}& \multirow{2}{*}{$\#$(patches)} & \multicolumn{3}{c|}{Polygonization}& \multicolumn{3}{c|}{Frame rate(fps)}& \multicolumn{2}{c}{Speedup} \\
  %\cline{5-10}\\
  &&&&sub&acc(ms)&lerp(ms)  &KSN08 &RS08 &Ours &$\frac{\mbox{Ours}}{\mbox{KSN08}}$ &$\frac{\mbox{Ours}}{\mbox{RS08}}$  \\

  \hline
  \ref{fig:drop}&$2\times2\times3$& $1\times1\times1$& 1&16         &0.67&0.28&11.59&13.10&57.62&4.97&4.39\\
  \ref{fig:stamp}&$2\times2\times2$& $1\times1\times3$&  3& 4       &0.94&0.54&23.62&27.44&97.65&4.13&3.55\\%78.12
  \ref{fig:01}&$2\times2\times2$& $17\times17\times17$  & 1636  &4  &6.96&1.82&10.66&11.54&49.57&4.65&4.29\\%19.83
  \ref{fig:04}&$2\times2\times2$& $13\times9\times9$& 566&4         &3.20&0.79&13.80&14.71&55.32&4.00&3.76\\%22.13
  \ref{fig:05}&$2\times2\times2$& $16\times16\times16$  & 2188  &8  &7.95&2.20&7.88 & 8.40&42.07&5.33&5.00\\%16.83
  \ref{fig:torus}&$2\times2\times2$& $13\times13\times5$& 248&4     &0.98&0.80&14.03&14.95&60.05&4.27&4.01\\%24.02
  \ref{fig:bishop}&$2\times2\times2$& $14\times24\times24$&1185& 4  &6.88&2.32&28.26&32.19&86.17&3.04&2.67\\%34.47
  \ref{fig:333}&$3\times3\times3$& $4\times18\times18$& 648&4       &5.87&0.88&11.12&12.67&50.77&4.56&4.00\\%20.31
  \ref{fig:L01}&$2\times2\times2$& $19\times15\times15$  & 1145& 4  &6.39&1.71&18.13&19.67&70.50&3.88&3.58\\%28.20
  \ref{fig:L04}&$2\times2\times2$& $9\times20\times20$& 1376&  1    &4.66&1.32&23.69&25.11&78.07&3.29&3.10\\%31.23
  \ref{fig:L02}&$2\times2\times2$& $37\times38\times38$& 3622& 1    &6.25&1.35&16.78&17.90&49.75&2.96&2.77\\%19.90
  \ref{fig:L05}&$2\times2\times2$& $23\times19\times19$& 1063& 1    &3.97&1.13&28.69&30.80&92.35&3.21&2.99\\%36.94
  \ref{fig:L03}&$2\times2\times2$& $18\times18\times18$& 937&  4    &7.52&2.27&16.66&18.20&66.47&3.98&3.65\\%26.59
  \ref{fig:L06}&$2\times2\times2$& $19\times19\times19$& 1494& 4    &7.27&1.81&13.13&14.08&56.17&4.27&3.98\\%22.47
  \hline
\end{tabular}
\end{table*}



\begin{figure}[h]
\centering
\includegraphics[width=0.45\textwidth]{Iter}
\caption{The numbers of function evaluations in our method on the
left and in the \textbf{RS08} method on the right, while the
\textbf{RS08} and \textbf{KSN08} methods have similar convergence
properties.} \label{fig:iterNum}
\end{figure}
We implemented the \textbf{KSN08} \cite{Kanamori08} method and the
\textbf{RS08} \cite{Reimers08} method as a comparison. The rendering
times and the speedup of our methods are shown is Table
\ref{table:patch}. The \textbf{RS08} method is extended to render
the ABS surface by adopting depth peeling as in \textbf{KSN08}.
According to our experiments, the B{\'e}zier clipping
method\cite{Kanamori08} is slightly slower than the knot insertion
method\cite{Reimers08} due to an increase in the number of
operations in each step. Aside from root-finding, the functional
composition itself will consume many function evaluations of the ABS
surface $F$, the surface normal $\nabla F$ and the univariate
polynomial. These are the most time-consuming parts in both our
method and functional composition methods. The color maps in
Fig.\ref{fig:iterNum} show the number of these function evaluations
in their methods and ours before a hitting point is determined for
an ABS surface of degree $6$. Their methods involve many more
function evaluations than ours. As a result, our algorithm achieves
a 3-5-fold speedup as shown in Table \ref{table:patch}. Even when
runtime of the isotopic polygonization in the preprocess step is
taken into account, our method is faster than the functional
composition methods. Two key reasons for the speedup are the initial
guesses provided by the rasterization of the isosurface and the root
isolation by adopting polar surfaces, and both of these can save us
from heavy overhead functional composition for piecewise algebraic
surfaces, because the functional composition contains large number
of the function evaluations.

% We can draw that the approaches will introduce heavy overhead for piecewise
%algebraic surface patches rendering. The speedup owes to  the
%initial guesses provided by rasterization of isosurface . In this
%way, we can purge away the non-intersection domains efficiently.

%For ray casting ABS surface, the main difference between the
%\textbf{KSN08} method and the \textbf{RS08} method are their
%root-finding methods of univariate polynomial. In the \textbf{RS08}
%method, the root is numerically solved by iteratively knot
%insertion~\cite{Mrken07}. The root-finding method in the
%\textbf{KSN08} method is B{\'e}zier clipping.


\section{Conclusion and Future Work}\label{sec:Conclusion}
In this paper, we presented a GPU-based real-time ABS surface
rendering algorithm, which performs ray and surface intersection
computations in object space. According to regular termination
criteria of domain subdivision of an ABS surface, an enhanced
polygonization is performed on GPU, where the isosurface has
accurate vertices and is isotopic to the ABS surface. The
rasterization result of the isosurface via traditional rendering
pipelines provides a rough initial guess for numerically solving for
the ray and ABS surface intersection. The silhouette, defined as the
intersection between the ABS surface and its polar surface is
scan-converted in floating-point accuracy via the Newton-Raphson
method. Then, the initial guess in the gap region between the ABS
surface and isosurface can be obtained via linear interpolation. The
polar surface is also employed to isolate the root near the
silhouette so as to improve the robustness of the root-finding
procedure. Apart from the silhouette and its vicinity, most of the
hit points between the ray and ABS surface can be solved via the
quadratically convergent Newton-Raphson method. Furthermore, the
surface normal that is important during illumination computation is
the byproduct of the Newton-Raphson method. We also design an
efficient and flexible data access to facilitate the performance on
GPU via the programming language CUDA. The experimental results show
that the proposed method can render the ABS surfaces of degree 6-9
composed of thousands of patches in real time.

\begin{figure}[h]
%\centering
\includegraphics[width=0.43\textwidth]{Error}
\caption{Artifacts caused by ill-conditioned
computations}\label{fig:Error}
\end{figure}
Several aspects of the proposed method can be improved in future
research. In the isotopic polygonization, each patch corresponds to
a grid resolution. We should choose the maximum resolution so that
the isosurface is free of cracks. This will lead to additional
costs. Aliasing artifacts may arise due to ill-conditioned
calculations. For example, when one patch is almost planar, and is
viewed from an extremely flat angle, silhouette computation will be
very sensitive to numerical noise as shown in Fig.\ref{fig:Error}.
Currently, other root-finding methods also struggle to handle this
case well because the nonlinear system will be under-determined.
However, the region with artifacts is very small in general.
Nevertheless, we are looking forward to seeking a more robust and
precise root-finding algorithm to address this problem. The criteria
used to determine whether a polar surface is null or not is
conservative. As a result, some overheads for root isolation are
required even in the case of a single root. The situation can be
improved by use of a more accurate function estimation.



\textbf{Acknowledgment} The research is funded by the NSF of China (60933007,
60873046), the 973 program of China (2009CB320801).
We would like to thank Dr. Aaron Knoll for
his valuable suggestions, which clarify the advantages and
difficulties in root-finding.

%Demo video can be downloaded at:
%\url{http://www.cad.zju.edu.cn/home/jqfeng/RealtimeRaycastingABS.htm}


%\newpage

\begin{figure*}[htb]
    \centering
    %\setcounter{subfigure}{2}
    \subfigure[1636 patches, $C^1$, 49.57 fps]
    {
        \includegraphics[width=0.4\textwidth]
        {01}
        \label{fig:01}
    }
    \subfigure[566 patches, $C^1$, 55.32 fps]
    {
        \includegraphics[width=0.4\textwidth]
        {04}
        \label{fig:04}
    }
    \subfigure[2188 patches, $C^1$, 42.07 fps]
    {
        \includegraphics[width=0.4\textwidth]
        {05}
        \label{fig:05}
    }
    \subfigure[248 patches, $C^1$, 60.05 fps]
    {
        \includegraphics[width=0.4\textwidth]
        {torus}
        \label{fig:torus}
    }
    \subfigure[1185 patches, $C^1$, 86.17 fps]
    {
        \includegraphics[width=0.4\textwidth]
        {bishop}
        \label{fig:bishop}
    }
    \subfigure[648 patches, $C^2$, 50.77 fps]
    {
        \includegraphics[width=0.4\textwidth]
        {333}
        \label{fig:333}
    }
    \caption{Rendering results. \subref{fig:333} is globally $C^{2}$ smooth, and all others are $C^{1}$ smooth.}
    \label{fig:demo}
\end{figure*}
\begin{figure*}[ht]
   \centering
    %\setcounter{subfigure}{8}
    \subfigure[1145 patches, $C^1$, 70.50 fps]
    {
        \includegraphics[width=0.35\textwidth]
        {01_15}\label{fig:L01}
    }
    \subfigure[1376 patches, $C^1$, 78.07 fps]
    {
        \includegraphics[width=0.45\textwidth]
        {01_03}\label{fig:L04}
    }
    \subfigure[3622 patches, $C^1$, 49.75 fps]
    {
        \includegraphics[width=0.35\textwidth]
        {01_02}\label{fig:L02}
    }
    \subfigure[1063 patches, $C^1$, 92.35 fps]
    {
        \includegraphics[width=0.45\textwidth]
        {01_04}\label{fig:L05}
    }
    \subfigure[937 patches, $C^1$, 66.47 fps]
    {
        \includegraphics[width=0.35\textwidth]
        {01_09}\label{fig:L03}
    }
    \subfigure[1494 patches, $C^1$, 56.17 fps]
    {
        \includegraphics[width=0.45\textwidth]
        {01_05}\label{fig:L06}
    }
    \caption{Rendering results. The models are constructed by approximation of triangle meshes and are all $C^{1}$ smooth.}
    \label{fig:demo3}
\end{figure*}
\appendix\section{Proof of Root isolation via polar surface}
Let ${\mathbf E}(x_e, y_e, z_e)$ and $\mathbf{d}(d_{x},d_{y},d_{z})$
be a viewpoint and ray direction respectively. Then, the ray
originated from $\mathbf E$ can be expressed as:
\begin{equation}\label{equ:ray}
\left \{
\begin{array}{ccc}
% \nonumber to remove numbering (before each equation)
  x(t) & = & x_e+t\cdot d_x  \\
  y(t) & = & y_e+t\cdot d_y  \\
  z(t) & = & z_e+t\cdot d_z
\end{array}
\right .
\end{equation}
The ray and surface intersection can be obtained by solving the
following univariate polynomial equation:
\begin{equation}\label{equ:univariate}
f(t)=F(x(t),y(t),z(t))=0.
\end{equation}
Its derivative polynomial is
\begin{eqnarray}\label{equ:derivative}
f'(t)&=&F_x\cdot d_x+F_y\cdot d_y+F_z\cdot d_z  \nonumber \\
     &=&F_x(\displaystyle\frac{x-x_e}{t})+F_y(\displaystyle\frac{y-y_e}{t})+F_z(\displaystyle\frac{z-z_e}{t}) \nonumber \\
     &=&\displaystyle\frac{nF-w_eF_w-x_eF_x-y_eF_y-z_eF_z}{t}\\
     &=&\displaystyle\frac{nF-P}{t} \nonumber
\end{eqnarray}
where the third step is obtained according to Euler's law, which is
the relationship between a homogeneous polynomial of degree $n$ and
its derivative:
\begin{equation}\label{equ:Poly&Deriv}
F(x,y,z,w)=\frac{xF_x+yF_y+zF_z+wF_w}{n}.
\end{equation}
By using Equation~(\ref{equ:derivative}), the functional composition
between the ray~(\ref{equ:ray}) and the polar
surface~(\ref{equ:Polar}) can be written as:
\begin{equation}\label{equ:polaruni}
p(t) = nf(t)-tf'(t)
\end{equation}
This is a linear combination of $f(t)$ and its derivative $f'(t)$.

\begin{mytheorem}\label{theorem:1}
Assuming that a polynomial $f(t)$ has two adjacent single roots
$t_0$ and $t_1$ in an interval $[a, b]$. we define a new polynomial
as:
\begin{equation}\label{equ:gx}
g(t)=c_0f(t)+c_1f'(t),
\end{equation}
where $c_0$ and $c_1$ are non-zero real numbers. Then for any given
$c_0$ and $c_1$, $g(t)$ has at least one root in the interval $(t_0,
t_1)$.
\end{mytheorem}

\begin{proof}
Let the root of $g(t)$ in the interval $(t_0, t_1)$ be noted as
$r(c_0, c_1)$. When $c_1\rightarrow 0$, $r(c_0, c_1)$ can approach
its lower and upper bounds $t_0$ and $t_1$ respectively.

If $t_0\neq t_1$, then according to Rolle's theorem, $f'(t)$ has at
least one root $t_2\in[t_0, t_1]$. Without loss of generality, we
assume $f(t_2)>0$. There exist $t_3\in[t_0, t_2]$ and $t_4\in[t_2,
t_1]$ such that $f(t_3)=f'(t_3)$ and $f(t_4)=-f'(t_4)$. Let
$k=c_1/c_0$, because $f(t)$ and $f'(t)$ are continuous, there is at
least one root $r(c_0, c_1)$ for $g(t)=0$ in
\begin{enumerate}
    \item   $[t0, t3]$ if $k\in(0, -1]$;
    \item   $[t3, t2]$ if $k\in[-1, -\infty)$;
    \item   $[t2, t4]$ if $k\in(\infty, 1]$;
    \item   $[t4, t1]$ if $k\in[1, 0)$;
\end{enumerate}
It is similar to the case of $f(t_2)>0$. Therefore, $g(t)$ has at
least one root in the interval $(t_0, t_1)$.\end{proof}

%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix

%% \section{}
%% \label{}

%% References
%%
%% Following citation commands can be used in the body text:
%% Usage of \cite is as follows:
%%   \cite{key}          ==>>  [#]
%%   \cite[chap. 2]{key} ==>>  [#, chap. 2]
%%   \citet{key}         ==>>  Author [#]

%% References with bibTeX database:

\bibliographystyle{model1b-num-names}
%\bibliographystyle{model1a-num-names}%use this template will cause error!
\bibliography{ABS_CAG}

%% Authors are advised to submit their bibtex database files. They are
%% requested to list a bibtex style file in the manuscript if they do
%% not want to use model1a-num-names.bst.

%% References without bibTeX database:

% \begin{thebibliography}{00}

%% \bibitem must have the following form:
%%   \bibitem{key}...
%%

% \bibitem{}

% \end{thebibliography}



\end{document}

%%
%% End of file `elsarticle-template-1a-num.tex'.
